close all
clear
clc

% problem: design a anti-reflection layers for 2.4GHz,1ft thick, ER=12 slab
% to make it, add ER2 material at both side of the slab
% so material structure would be |ER2|ER_slab|ER2|

% Common value for simulation
c0 = 299792458; % light speed
gigahertz = 1e9;
ER_src = 1; % ER at source
UR_src = 1; % UR at source
ER_air = 1;
UR_air = 1;
nmax = sqrt(ER_air*UR_air); % refractive index for material
nbc = 1; % refractive index at the grid boundaries
nsrc = 1; % refractive index at the soruce cell

freq_max = 1*gigahertz; % 1GHz

x_length = 1; % 1m in x axis
y_length = 2; % 1m in y axis


% =================================================================
% Start calculate grid resolution
NRES = 20; % resolve the wave with at least 10 cells, better >=10
lambda_min = c0./freq_max./nmax;
delta_lambda = lambda_min/NRES;

% use free space in this example so material resultion is not needed
%NDRES = 4; % normally 1~4 for resolution of feature size
%delta_d = x_length/NDRES;  
%d_temp = min(delta_lambda, delta_d);

d_temp = delta_lambda;

Nx = ceil(x_length/d_temp);
Ny = ceil(y_length/d_temp);
dx = x_length/Nx; % get the size of grid
dy = y_length/Ny; % get the size of grid
d_avg = (dx+dy)/2;


% =================================================================
% Start calculate delta_t and tau 
dt = nbc*d_avg/(2*c0);
tau = 1/(2*freq_max);
t0 = 6*tau;
t_prop = max(nmax*Nx*dx/c0, nmax*Ny*dy/c0);

T_total = 12*tau + 5*t_prop;
STEPS = ceil(T_total/dt);
% =================================================================

ER = ones(Nx,Ny);
UR = ones(Nx,Ny);

% =================================================================
% Start calculate Source for Ez mode
t=(0:STEPS-1).*dt; % time axis
delta_t = nsrc*d_avg/(2*c0) + dt/2; % total delay between E and H
A = -sqrt(ER_src/UR_src); % amplitude of H field
Esrc = exp(-((t-t0)/tau).^2); % E filed source
Hsrc = A*exp(-((t-t0+delta_t)/tau).^2); % H field source
% =================================================================

% Compute updated coefficients
mUR_xx = -(c0*dt./UR);
mUR_yy = -(c0*dt./UR);
mHz = dt*c0;

% Initialize field to zero
Ez = zeros(Nx,Ny); 
CEx = zeros(Nx,Ny); 
CEy = zeros(Nx,Ny); 
CHz = zeros(Nx,Ny); 
Hx = zeros(Nx,Ny); 
Hy = zeros(Nx,Ny); 
Dz = ER.*Ez;

fig=figure;
set(fig,'Name', 'FDTD 2D Free Space Simulation');
set(fig,'NumberTitle', 'off');

% position of the source, 2nd cell
Nx_src = 20;
Ny_src = 30;

xa=1:Nx;
ya=1:Ny;
[Y,X] = meshgrid(ya,xa);

% Main FDTD Loop for Ez mode
for T = 1 : STEPS
    
   % Update CEx
   for nx = 1 : Nx
       for ny = 1 : Ny-1
           CEx(nx,ny) = (Ez(nx,ny+1) - Ez(nx,ny))/dy;
       end
       CEx(nx,Ny) = (0 - Ez(nx,Ny))/dy; % Dirichlet Boundary Conditions
   end
   % Update CEy
   for ny = 1 : Ny
       for nx = 1 : Nx-1
           CEy(nx,ny) = -(Ez(nx+1,ny) - Ez(nx,ny))/dx;
       end
       CEy(Nx,ny) = -(0 - Ez(Nx,ny))/dx; % Dirichlet Boundary Conditions
   end
   % Update H from E
   Hx = Hx + mUR_xx .* CEx;
   Hy = Hy + mUR_yy .* CEy;

   
   % Update CHz
   CHz(1,1) = (Hy(1,1) - 0)/dx...
             -(Hx(1,1) - 0)/dy;
   for nx = 2:Nx
      CHz(nx,1) = (Hy(nx,1) - Hy(nx-1,1))/dx...
                 -(Hx(nx,1) - 0)/dy;
   end
   for ny = 2:Ny
       CHz(1,ny) = (Hy(1,ny) - 0)/dx...
                  -(Hx(1,ny) - Hx(1, ny-1))/dy;
       for nx = 2:Nx
           CHz(nx,ny) = (Hy(nx,ny) - Hy(nx-1,ny))/dx...
                       -(Hx(nx,ny) - Hx(nx, ny-1))/dy;
       end
   end
     
   % update D from H
   Dz = Dz + mHz .* CHz;
   % update source
   Dz(Nx_src,Ny_src) = Dz(Nx_src,Ny_src) + Esrc(T)*ER(Nx_src,Ny_src);
   % update E from D
   Ez = Dz./ER;
   
   %Plot
   pcolor(X,Y,abs(Ez).^2);
   axis ij
   axis equal tight;
   colormap('jet');
   shading interp;
   %caxis([0,1]);
   %colorbar
   hold on
   title(sprintf('Step: %d of %d',T, STEPS));

   drawnow;
   hold off
   %pause(0.01);
     
end









